
<?php
function ecuts() {
        $ev = array(1e-100,1e-10,1e-5,0.001,0.1,1,10,100);
        return $ev;
}

// remove non-alphabetical characters
function blastseqpack($seq) {
    $seq = trim($seq);
    if (strlen($seq)==0) return $seq;
    
    $lines = explode("\n",$seq);
    $line = trim(current($lines));
    while (empty($line)) {
            next($lines);
            $line = trim(current($lines));
    }
    if (substr($line,0,1)==">") {
        next($lines);
        $line =trim(current($lines));
    }
    
    do {
        //$return_this = $return_this.preg_replace("/^[a-zA-Z]/","",$line);
        $return_this = $return_this.eregi_replace("[^a-z]","",$line);
        next($lines);
        $line = trim(current($lines));
    } while ($line && !empty($line) && substr($line,0,1)!=">");

    $return_this = strtoupper($return_this);
    return $return_this;
}

// parse blast result
function blastparse($blastout) {
    $lines = file($blastout);
    while (current($lines)) {            
        next($lines);
        $line = trim(current($lines));
        if (strstr(line,"***** No hits found ******")) {
            return null;
        }
        if (strstr($line,"Sequences producing significant alignments:")) {
            next($lines);
            next($lines);
            while (current($lines)) { // skip hit list
                if (trim(current($lines))=="") {
                    break;
                }
                next($lines);
            }
            break;
        }
    }

    // hit data begins
    $hits = array();
    while (current($lines)) {
        while (true) {
            $line = current($lines);
            if (!$line) { // end of file
//print_r($hits);
                return $hits;
            }
            if (trim($line)!="") {
                if (strstr($line,">")) {
                    break;
                }
            }
            next($lines);
        }
        // one hit
        $hit = array();
        $line = trim(current($lines)); //>dds#pps|1|10
        $strs = explode('|', substr($line,1));
        $hit["nrseq"] = $strs[0];
        $hit["protein"] = $strs[1];
        $hit["location"] = $strs[2];
        
       // while (strpos(current($lines)," Score = ")==0) {
        //    next($lines);
       // }
        next($lines);
        next($lines);
        next($lines);
        
        $line = trim(current($lines)); // Score = 54.3 bits (129), Expect = 3e-010
        $s = strpos($line,"Score = ");
        $line = trim(substr($line,$s+8));
        $s = strpos($line," bits (");
        $hit["score"] = substr($line, 0, $s);
        $line = trim(substr($line,$s+7));
        $s = strpos($line,"), Expect = ");
        $hit["bits"] = trim(substr($line, 0, $s));
        $hit["evalue"] = floatval(substr($line,$s+12));
        
        next($lines);
        $line = trim(current($lines)); // Identities = 10/22 (45%), Positives = 13/22 (59%), Gaps = 3/22 (13%)
        $arr = explode(")",$line);
        $hit["ide"] = substr($arr[0],strpos($arr[0],"(")+1);
        $hit["pos"] = substr($arr[1],strpos($arr[1],"(")+1);
        if (count($arr==2)) {
            $hit["gap"] = "0%";
        } else {
            $hit["gap"] = substr($arr[2],strpos($arr[2],"(")+1);
        }
        
        next($lines);
        next($lines); 
        
        $line = trim(current($lines)); //Query: 184 HSLHSYFLLVGD 195
        $s = strpos($line,"Query:");
        $line = trim(substr($line,$s+6));
        $s = strpos($line," ");
        $hit["querystart"] = intval(substr($line,0,$s))-1;
        $line = trim(substr($line, $s));
        $s = strrpos($line," ");
        $hit["queryseq"] = trim(substr($line,0,$s));
        $hit["queryend"] = intval(trim(substr($line,$s)))-1;
        
        next($lines); 
        $line = trim(current($lines)); //T  D+P+FG+  G   GQ  GQ
        $hit["matchseq"] = $line;
        
        next($lines);
        $line = trim(current($lines));
        $s = strpos($line,"Sbjct:");
        $line = trim(substr($line,$s+6));
        $s = strpos($line," ");
        $hit["sbjctstart"] = intval(substr($line,0,$s))-1;
        $line = trim(substr($line, $s));
        $s = strrpos($line," ");
        $hit["sbjctseq"] = trim(substr($line,0,$s));
        $hit["sbjctend"] = intval(trim(substr($line,$s)))-1;

        // index of site are ungapped
        $hit["sites"] = getsbjctsites($hit["nrseq"], $hit["protein"], $hit["location"], $hit["sbjctstart"], $hit["sbjctend"]);

        if (count($hit["sites"])>0) {
            $hits[] = $hit;            
        }

        next($lines);
    }
}

function getsbjctsites($nrseq, $protein, $location, $sbjctstart, $sbjctend) {
        $locInPep = getSiteFromNrPepSeq($nrseq);
        $sites = array();
        foreach ($locInPep as $loc) {
            if($loc<$sbjctstart||$loc>=$sbjctend)
                continue;
            $siteLoc = $loc + $location;
            $query = "SELECT id FROM site WHERE protein=$protein AND locationInProtein=$siteLoc AND ptm='p'";
            $result = mysql_query($query) or die("Error, query failed: $query");
            while($row = mysql_fetch_array($result, MYSQL_ASSOC)) {
                $sites[$loc-$sbjctstart] = $row['id'];
            }
        }
        return $sites;
}

function getSiteFromNrPepSeq($nrpseq) {
    $curr = strpos($nrpseq, '#');
    while ($curr) {
        $sites[] = $curr - count($sites) - 1;
        $curr = strpos($nrpseq, '#', $curr+1);
    }
    return $sites;
}

//// get sites for hit
//function getsbjctsites($id, $sbjctstart, $sbjctend) {
//        $sites = array();
//        $query = "SELECT siteid,loc FROM peptidesiterelation WHERE pepid = $id ORDER BY loc";
//        $result = mysql_query($query) or die("Error, query failed: $query");
//        while($row = mysql_fetch_array($result, MYSQL_ASSOC)) {
//                $loc = $row["loc"];
//                if ($loc<$sbjctstart||$loc>=$sbjctend) {
//                        continue;
//                }
//                $index = $loc-$sbjctstart;
//                $siteid = $row["siteid"];
//                $sites[$index] = $siteid;
//        }
//        //print_r($sites);
//        return $sites;
//}

// align blast result
function blastalign(&$query,&$hits) {

    // recode the gap "-" first
    $insert_query = array(); // gap inserted to query
    $inserts_sbjct = array(); // gap inserted to subjects
    foreach ($hits as $i=>$hit) {
        $queryhit = $hit["queryseq"];
        $querystart = $hit["querystart"];

        $insert_sbjct = array(); // gap inserted to subject
        $index = strpos($queryhit,"-");
        while ($index!==false) {
            $insert_sbjct[$index]++;
            $queryhit =  substr($queryhit,0,$index).substr($queryhit,$index+1); //remove '-'
            $index = strpos($queryhit,"-");
        }

        if (count($insert_sbjct)>0) {
            $inserts_sbjct[$i] = $insert_sbjct; //keep for use later
        }

        // for the same position, put the maximum number of gaps
        foreach ($insert_sbjct as $index=>$num) {
            $index += $querystart;
            if ($num>$insert_query[$index]) {
                $insert_query[$index] = $num;
            }
        }
    }
//print_r($insert_query);
//print_r($inserts_sbjct);

    // insert "-" into the query sequence
    if (count($insert_query)>0) {
        krsort($insert_query,SORT_NUMERIC);
        foreach($insert_query as $index=>$num_query) {
            $insert = str_repeat("-",$num_query);
            $query = substr($query,0,$index).$insert.substr($query,$index);
        }
    }

    foreach ($hits as $i=>$hit) {
        $queryseq = &$hit["queryseq"];
        $sbjctseq = &$hit["sbjctseq"];
        $matchseq = &$hit["matchseq"];

        $querystart = $hit["querystart"];
        $queryend = $hit["queryend"];
        $insert_sbjct = $inserts_sbjct[$i];

        //for ($index=$querystart; $index<$queryend; $index++) {
        foreach ($insert_query as $index=>$num_query) {
                //$num_query = $insert_query[$index];
                if ($querystart<=$index && $index<=$queryend) {
                        $index -= $querystart;
                        $num_sbjct = $insert_sbjct[$index];
                        if ($num_query>$num_sbjct) { // more gaps need to be inserted
                            $insert = str_repeat("-",$num_query-$num_sbjct);
                            $index_gapped = indexConvert_gapExcluded($queryseq,$index);
                            if ($index_gapped!=0) {
                                $queryseq = substr($queryseq,0,$index_gapped).$insert.substr($queryseq,$index_gapped);
                                $sbjctseq = substr($sbjctseq,0,$index_gapped).$insert.substr($sbjctseq,$index_gapped);
                                $matchseq = substr($matchseq,0,$index_gapped).$insert.substr($matchseq,$index_gapped);
                            }
                        }
                }
        }

        $querystart = strpos($query,$queryseq);
        if ($querystart===false) {
            echo "error:$queryseq<br>";
            echo $query;
            $hit["error"]==true;
            print_r($hit);
        }

        $queryend = $querystart+strlen($queryseq);

        $hit["querystart"] = $querystart;
        $hit["queryend"] = $queryend;

        $hits[$i] = $hit;
        // note : $hit["sites"] has not been syncronised here, caution when use
    }

}

function group_hit($query,$hits) { // no use
    $hits_group = array();// array of array of hit

    // find a position for each hit
    foreach($hits as $hit) {
        $start = $hit["querystart"];
        $end = $hit["queryend"];

        foreach($hits_group as $i=>$hits_it) {
            $conflict = false;
            foreach($hits_it as $hit_it) {
                if (!($end<$hit_it["querystart"]||$start>$hit_it["queryend"])) { // conflict
                    $conflict = true;
                    break;
                }
            }

            if (!$conflict) {
                $hits_it[$start] = $hit;
                $hits_group[$i] = $hits_it;
                break;
            }

            $hits_group[$i] = $hits_it;
        }

        // if no suitable position
        if ($conflict||empty($hits_group)) {
            $hits_new[$start] = $hit;
            $hits_group[] = $hits_new;
        }
    }

    foreach($hits_group as $i=>$hits_it) {
        ksort($hits_it,SORT_NUMERIC); // sort asc
        $hits_group[$i] = $hits_it;
    }

    return $hits_group;
}


function blastformat($query,$hits,$evalues,$aaperline,$org,$ref) {
//print_r($hits_group);

   // format
   $naa = strlen($query);
   $nline = ceil($naa/$aaperline);
   $seqarrs = array();

//print_r($seqarrs);

   //$nhit = 0;

   $querysites = array();
   // subject
   foreach($hits as $hit) {
           $nrseq = $hit["nrseq"];
           $protein = $hit["protein"];
           $location = $hit["location"];
           $evalue = $hit["evalue"];
           $color = getcolor($evalue,$evalues);

           $start = $hit["querystart"];
           $end = $hit["queryend"];
           
           $sbjct = $hit["sbjctseq"];
           $sbjctstart = $hit["sbjctstart"];
           $sbjctend = $hit["sbjctend"];
           $sbjctsites = $hit["sites"];

           //$sites = getsbjctsites($id, $sbjct, $sbjctstart, $sbjctend);
           //if (count($sites)==0) continue; //skip if site not in

           //$nhit++;

           //record for query
           setquerysites($querysites,$sbjct,$sbjctsites,$start,$org,$ref);
           
           $il = floor($start/$aaperline);
           $s = $il*$aaperline;
           while ($end>$s) {
               $t = min($s+$aaperline,$naa);
               $ss = max($start,$s);
               $tt = min($end,$t);
               $seq = str_repeat("&nbsp;",$ss-$s); // fill the blank
               
               if ($ss==$start) $seq = $seq."&lt;";
               else $seq = $seq.$seq."&nbsp;";


               $prea = formata($nrseq,$protein,$location,$evalue,$color, $sbjctsites,$org,$ref);
               $seq = $seq.$prea;
               
               //$seq = $seq.substr($sbjct,$ss-$start,$tt-$start);
               $seq = $seq.formatpepsbjct($sbjct,$sbjctsites,$ss-$start,$tt-$start,$prea,$org,$ref);
               $seq = $seq."</a>";

               if ($tt==$end) $seq = $seq."&gt;";

               $seqarrs[$il]["sbjct"][] = $seq; // add

               $il++;
               $s = $il*$aaperline;
           }
   }

        ksort($querysites);
        //print_r($querysites);
   // query
   for ($il=0; $il<$nline; $il++) {
       $s = $il*$aaperline;
       $t = min($s+$aaperline,$naa);

       //$seqarrs[$il]["query"] = substr($query,$s,$t-$s);
       $seqarrs[$il]["query"] = formatpepquery($query,$querysites,$s,$t);
   }

//print_r($seqarrs);
   ksort($seqarrs);

   return $seqarrs;

   //$ret["nhit"] = $nhit;
   //$ret["seqarrs"] = $seqarrs;
   //return $ret;

}

function formatpepquery($seq,$sites,$start,$end,$org,$ref) {
        $ret = "";
        foreach ($sites as $index=>$siteids) {
                if ($index<$start) continue;
                if ($index>=$end) break;

                $ps = substr($seq, $index, 1);
                if (!preg_match("/[stySTY]/",$ps)) continue;

                $ret = $ret.substr($seq, $start, $index-$start);

                //print_r($siteids);
                $siteids = array_unique($siteids);
                $ids = implode(",",$siteids);
                $highlight = "highlight('site".implode("');highlight('site",$siteids)."');";
                $dehighlight = "dehighlight('site".implode("');dehighlight('site",$siteids)."');";
                $ret = $ret.
//                       "<a href=\"../data/phosphosite.php?id=$ids\"".
                       "<a href=\"#\"".
                       " style=\"color:blue\"".
                       " onmouseover=\"$highlight\"".
                       " onmousedown=\"$dehighlight\"".
                       " onmouseout=\"$dehighlight\"".
                       ">$ps</a>";
                $start = $index+1;
        }
        $ret = $ret.substr($seq,$start,$end-$start);
        return $ret;
}

function formata($nrseq,$protein,$location,$evalue,$color,$sbjctsites,$org,$ref) {
        $pep_node_id = "pep.$protein.$location.".str_replace('#', '_', $nrseq);
        $highlight = "highlight('$pep_node_id');";
        $dehighlight = "dehighlight('$pep_node_id');";
        foreach ($sbjctsites as $index=>$siteid) {
                $highlight = $highlight."highlight('site$siteid');";
                $dehighlight = $dehighlight."dehighlight('site$siteid');";
        }

        $format_html = str_replace('#', '%23', $nrseq);
        $a = "<a href=\"../data/peptide.php?pnrseq=$format_html&pro=$protein&loc=$location&org=$org&ref=$ref\"".
                      " title=\"E-value:$evalue\"".
                      " name=\"pep$id\"".
                      " onmouseover=\"$highlight\"".
                      " onmousedown=\"$dehighlight\"".
                      " onmouseout=\"$dehighlight\"".
                      " style=\"color:$color\"".
                      ">";

        return $a;
}

//
function formatpepsbjct($seq,$sites,$start,$end,$prea,$org,$ref) {
        $ret = "";
        foreach ($sites as $index=>$siteid) {
                $index = indexConvert_gapExcluded($seq, $index);
                if ($index<$start) continue;
                if ($index>=$end) break;
                $ret = $ret.substr($seq, $start, $index-$start);
                $ret = $ret."</a>".
                        "<a name=\"site$siteid\"".
                        "   href=\"../data/phosphosite.php?id=$siteid&org=$org&ref=$ref\"".
                        "   onmouseover=\"highlight('site$siteid')\"".
                        "   onmousedown=\"dehighlight('site$siteid')\"".
                        "   onmouseout=\"dehighlight('site$siteid')\"".
                        "   style=\"color:blue\">".
                        substr($seq, $index, 1).
                        "</a>$prea";
                $start = $index+1;
        }
        $ret = $ret.substr($seq,$start,$end-$start);
        return $ret;
}

//// get sites for hit
//function getsbjctsites($id, $sbjct, $sbjctstart, $sbjctend) {
//        $sites = array();
//        $query = "SELECT siteid,loc FROM peptidesiterelation WHERE pepid = $id ORDER BY loc";
//        $result = mysql_query($query) or die('Error, query failed');
//        while($row = mysql_fetch_array($result, MYSQL_ASSOC)) {
//                $loc = $row["loc"];
//                if ($loc<$sbjctstart||$loc>=$sbjctend) continue;
//                $index = $loc-$sbjctstart;
//                $index = indexConvert_gapExcluded($sbjct, $index);
//                $siteid = $row["siteid"];
//                $sites[$index] = $siteid;
//        }
//        return $sites;
//}

function setquerysites(&$querysites,$sbjctseq,$sbjctsites,$querystart) {
        foreach($sbjctsites as $index=>$siteid) {
                $index = indexConvert_gapExcluded($sbjctseq, $index);
                $querysites[$index+$querystart][] = $siteid;
        }
}


function getevalues($hits) {
        $evalues = array();
        foreach($hits as $hit) {
           $evalues[] = $hit["evalue"];
        }

        $evalues = array_unique($evalues);
        sort($evalues);

        return $evalues;
}

function getcolor($evalue,$evalues) {
        $n = count($evalues);

        $i = 0;
        foreach($evalues as $ev) {
                if ($evalue<$ev) {
                        $i++;
                }
        }

        $t = 230;
        $s = 50;
        $b = ($t-$s)/$n;

        $c = floor($s+$b*$i);

        //return "#".dechex($c)."0000";
        return "#00".dechex($c)."00";
        //return "#0000".dechex($c);
        //return "#".dechex($c)."00".dechex($c);
        //return "#00".dechex($c).dechex($c);

}

function color_legend($evalues) {
        //$evalues = array_reverse($evalues);
        $legend = "<table width=\"300\"><tr>";
        $evalue = current($evalues);
        $color = getcolor($evalue,$evalues);
        $legend = $legend."<td width=\"50\" bgcolor=\"$color\"><font color=\"white\" size=\"1\">".$evalue."</font></td>";
        foreach($evalues as $evalue) {
                $color = getcolor($evalue,$evalues);
                $legend = $legend."<td bgcolor=\"$color\"></td>";
        }
        $legend = $legend."<td align=\"right\" width=\"50\" bgcolor=\"$color\"><font color=\"white\" size=\"1\">".$evalue."</font></td>";
        $legend = $legend."</tr></table>";
        return $legend;
}


function indexConvert_gapExcluded($seq, $index) {
        if ($index<0 || $index>=strlen($seq)) return -1;

        $count = 0;
        for ($i=0; $count<$index; $i++) {
            //if (substr($seq,$i,1)!="-"){
            if ($seq[$i]!='-') {
                $count++;
            }
        }

        return $i;
}

?>